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Stability of the branching structure of an RNA molecule is an important condition for its function. 
In this letter we show that the melting thermodynamics of RNA molecules is very sensitive to their 
branching geometry for the case of a molecule whose groundstate has the branching geometry of a 
Cayley Tree and whose pairing interactions are described by the Go model. Whereas RNA molecules 
with a linear geometry melt via a conventional continuous phase transition with classical exponents, 
molecules with a Cayley Tree geometry are found to have a free energy that seems smooth, at least 
within our precision. Yet, we show analytically that this free energy in fact has a mathematical 
singularity at the stability limit of the ordered structure. The correlation length appears to diverge 
on the high-temperature side of this singularity. 

PACS numbers: 87.15.Aa, 64.60.Fr, 87.15. Cc, 87.15. Nn 



A fundamental principle of statistical mechanics 
states that phase transitions are not possible for one- 
dimensional systems unless long-range interactions are 
present. It thus came as a surprise when Poland and 
Scheraga (PS) showed that an infinite, linear molecule 
composed of two flexible polymer strands bound to- 
gether by a local attractive interaction does undergo a 
true phase transition at the temperature where the two 
strands separate. The required long-range correlations 
are due to the fact that the partition function of a strand 
separation "bubble" has a power-law dependence on size. 
The mean bubble size, the correlation length, diverges at 
the critical point if the transition is continuous Q • This 
observation was particularly interesting because that sys- 
tem could be viewed as a simple model for the denatura- 
tion of double-stranded B-DNA molecules. 

The PS mechanism can be extended to the melting of 
- more complex — RNA molecules In a biologi- 
cal context, RNA molecules usually operate in a single- 
stranded mode. This single strand can however bend onto 
itself so the bases of the strand can self-pair into a pat- 
tern of bubbles and "stems" that can be displayed in the 
form of a tree-like planar graph, the " secondary struc- 
ture" [J] . The minimum-energy secondary structure of a 
functional RNA molecule plays an important role in its 
functioning, and can be predicted from the primary se- 
quence of nucleotides [B| . Melting of a minimum-energy 
secondary structure produces a "molten globule" state 
with the molecule fluctuating over a range of different sec- 
ondary structures Q . Importantly, in this molten globule 
state, most bases remain paired in contrast to the fully 
denatured state, which is favored at higher temperatures, 
with most of bases unpaired. In his pioneering paper of 
1968 0, de Gennes showed that the partition function 
G(L) of a large RNA molecule fluctuating over all possi- 
ble secondary structures with identical non-specific pair- 
ing energies has a power-law dependence on size of the 



FIG. 1: Single-stranded RNA molecule having a branched 
secondary structure that follows the outline of a Cayley Tree. 
Nucleotides are schematically indicated by circles, bonds be- 
tween nucleotides by a solid line and complementary pairing 
interactions by dashed lines, a) Groundstate structure with 
pairing restricted to a complementary "native" pair for each 
branch of the Cayley Tree, b) In a molten globule bubble 
(hatched) all possible pairing interactions are permitted. 



form ffr with 9 = 3/2. Bundschuh and Hwa |8j (BH) 
extended this result to show that if the groundstate of 
an RNA molecule is a long, linear hairpin stabilized by 
specific pairing energies then thermal fluctuations in the 
form of molten-globule bubbles produce a melting ther- 
modynamics that, formally, has the same form as that of 
the PS model. 

Actual RNA secondary structures have a branched, 
tree-like form, which raises the question how and if the 
melting thermodynamics of such a form differs from that 
of a simple hairpin. It is the experience with many sta- 
tistical mechanics models defined on tree-like geometries 
without circuits that they exhibit mean-field type critical 
behavior. Since, in the absence of excluded volume inter- 
actions, the critical properties of the PS and BH models 
already are of mean-field character, one would expect the 
free energy of branched secondary structures to exhibit 
mean-field critical behavior. In this letter we will show 
that in fact the melting thermodynamics of a particular, 
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highly branched secondary structure is highly anomalous: 
the numerically computed free energy appears to have no 
singularity, yet, surprisingly, we can demonstrate ana- 
lytically that the free energy does have a mathematical 
singularity at a point where the branched groundstate 
becomes unstable. The correlation length appears to di- 
verge on the high-temperature side of the singularity, yet, 
on the low-temperature side this singularity is not asso- 
ciated with a divergence of the correlation length. 

To demonstrate these claims, we consider an RNA 
molecule that has the shape of a Cayley Tree (see Fig.l). 
In the groundstate, the single strand traces out the 
perimeter of the tree, starting and ending at the root of 
the tree, with each branch of the tree occupied by a single 
complementary base-pair. The size of the molecule is in- 
dexed by the level k of the tree that is related to the total 
sequence length N(k) of the strand by N(k) = 2 k+2 - 2 
bases (a k = 1 tree is here a three-armed star with 
one base- pair per arm). After sequentially numbering 
the bases of the strand, one can denote this "designed" 
groundstate by a list S = {12,32}, ■ ■ ■ , {«m, 3m} 

of complementary pairs. We will assign a specific binding 
energy — e to any pair in this list. Pairing between two 
bases that do not appear in this list still will be allowed 
as long as it does not introduce any circuits (or "pseu- 
doknots") but the associated binding energy — e will be 
assumed to be less attractive than — e. This definition of 
the pairing energy, known as a "Go Model" guaran- 
tees that the secondary structure of the groundstate has 
the shape of a Cayley Tree. 

Our strategy to obtain the finite-temperature partition 
function of the system is to generalize the method of BH 
for the one-dimensional case by expressing the partition 
function in the form of a sum over all possible insertions 
of molten globule bubbles in the ground-state structure. 
Inserting a bubble into a Cayley Tree is more complex 
than into a linear structure: a bubble inside the tree can 
have different numbers of branches attached to it so one 
has to keep track of different bubble species. We will 
show elsewhere that, within the Go model, the partition 
function of any "designed" secondary structure can be 
written as a sum over configurations classified according 
to the size 2n of the "accessible" open bubble located at 
the base of the tree (see Fig.l). Here, n is the number of 
base pairs of the open bubble. Specifically, the partition 
function Z(k) can be written as: 

N(k)/2 

Z(k)= G(2n)W(k,n). (1) 

71=0 

In Eq. |T]), G(M) « is the partition function of a 
strand of length M with no specific pairing, i.e. all paired 
bases have a binding energy — e even if the pair appears in 
the list S of specific groundstate pairs. Next, W(k, n) is 
a restricted partition function, i.e., the partition function 
of a molecule with n accessible bases in the open bubble 



at the root, but not including the configurations of the 
open bubble. This restricted partition function can be 
written as a sum over all possible bubble insertions: 

W(k,n)= m lS ' in)l II (2) 

S'(n)cS {L(S')} 

Here, q — exp(f3e) and q — exp(/?e) while S'(n) is any of 
the subsets of S that is compatible with n base-pairs in 
the open bubble at the root. The number of specifically 
paired bases of S'(n) is denoted by |S"(n)|. Each term of 
Eq. @ represents a secondary structure having 
specifically paired bases linked together by a distribution 
of closed bubbles with sizes {L(S')}. 

Using Eq. ([2]) one can construct two linked recursion 
relations. First, cut a tree with restricted partition func- 
tion W(k,n) into two equal sized sub-trees with level 
index k — 1. The number of accessible base pairs of the 
two sub-trees together must add to n — 1 , as we removed 
one pair by the cutting operation. Because we permit no 
circuits, the restricted partition function of a level k tree 
and n > can be expressed in terms of a product of the 
restricted partition functions of two k — 1 level sub-trees: 

ra-1 

W(k,n) = ^2w(k-l,m)W(k-l,n-l-m) (3) 

m=0 

with W{k - l,m) = if m > 2 k - 1. The n = case 
— a tree with no bubble at the root — must be treated 
separately. Take the first complementary pair at the root 
of the tree out of the partition function, and then sum 
over all possible sizes for the bubble that immediately 
follows this pair (including a bubble of zero size). Now 
treat that bubble as the bubble at the root of a new tree 
that can again be cut into two equal parts in the same 
way as before. This leads to a second recursion relation: 

W(k,n)= (4) 

2 fc -l 2 fc -l 

(9-9E ^W(k~l,n 1 )W{k-l,n 2 )G{2(n 1 +n 2 )) 

ni—0 U2—0 

Equations (J3J) and ([!]) together constitute a complete 
set of recursion relations for W(k,n) that can be solved 
iteratively. The initial conditions for the recursion re- 
lations are W(1,0) = (q - q)[l + 4q + q 2 + 2q + q 2 }, 
W(l, 1) = (q- <z) 2 , W(l, 2) = 2{q- q), and W(l, 3) = 1, 
as follows by inspection. 

We carried out this iteration procedure numerically up 
to level k — 19 for different values of q = exp(/3e) and for 
fixed q = 4. In Fig. 2 we show the second derivative of the 
free energy per site with respect to q, which effectively 
correspond to the heat capacity. As one increases the 
value of k, a maximum develops near q. However, within 
the numerical precision, the free energy per site does not 
develop a thermodynamic singularity in the large N limit. 
This must be contrasted with the case where the molecule 
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has a linear hairpin groundstate, in which case the heat 
capacity very clearly develops such a singularity for much 
smaller system sizes (see inset of Fig. 2). 
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FIG. 2: Second derivative of the free energy with respect to 
the Boltzmann weight q of specifically paired bases plotted as 
a function of q for different values of the level k of the Cayley 
tree groundstate. The free energy was computed numerically 
from the recursion relations Eqs. © and Q and expressed 
in units of NksT with N the sequence length of the RNA 
strand. The arrow denotes the location of the mathematical 
singularity associated with melting of the root of the Cayley 
Tree. Inset: same except that the groundstate is a linear 
hairpin groundstate. A singularity develops near q = 18.4. 

In order to examine sub-leading contributions to the 
free energy, i.e., terms that are small compared to the 
leading term proportional to N, we also computed the 
"pinching free energy" 



AF(k) /k B T = In Z(k + 1) - 2 In Z(k) 



(5) 



For example, in a molten-globule phase the parti- 
tion function should have the asymptotic scaling form 
a + Zq /N 3 / 2 for large N. The pinching free energy 
AF(k)/k B T = §(fc + 2) In 2 - lna+ then would have a 
linear dependence on k, with slope 3/2. In an ordered 
phase, the partition function should scale as a~ Zq for 
large N, in which case AF(k)/ksT = — lna~ should be 
a constant independent of k. Figure 3 shows that, for q 
values up to 80, AF(k) indeed has a linear dependence 
on k for large k, with a slope close to 3/2 In 2. This indi- 
cates that, for q values below 80, the tree is in the molten- 
globule phase. Since for the corresponding case of a linear 
groundstate, the melting point is as low as q c = 18.4 for 
q = 4.0, we are forced to conclude that branching has a 
powerful destabilizing effect on the ordered state. 

For smaller k values, the pinching free energy is a con- 
stant, which indicates that the ordered groundstate dom- 
inates over shorter length scales. The crossover point be- 
tween the two regimes can be interpreted as a correlation 
length £ whose physical meaning would be that of the typ- 
ical size of smaller ordered Cayley Tree type structures 
imbedded in a larger molten-globule state. The value of 
£ increases with q according to Fig. 3 and beyond q = 80 
it exceeds our maximum system size (N = 10 6 ). A fit to 



a power-law £ ~ (q c — q) v produces a correlation length 
exponent v « 2.1 and a critical q c « 80. 




FIG. 3: Numerically computed "pinching" free energy AF(k) 
(see Eq. (J5J) versus the level k of the Cayley Tree. For q larger 
then 90, AF(k) /ksT is independent of k, consistent with the 
ordered groundstate. For q less than then 20, AF(k) /ksT can 
be fitted by the relation AF(k)/k B T « (fc + 2) ln2-lna+ for 
the molten globule state. The cross-over point between these 
two regimes for intermediate values of q marks the size of the 
ordered, correlated regions in the molten globule state. For q 
above 80, the size of the correlated regions exceeds the system 
size. 

Can we really be sure that there is a thermodynami- 
cally stable, ordered phase at low but finite temperatures 
or might the groundstate only appear at T = 0? In the 
ordered phase, the restricted partition function would be 
expected to scale as W(N, n) ~ w(h)zq asymptotically 
for large N. Here, w(n) is the fraction of configurations 
that have an open bubble at the root of size n. If we 
insert this Ansatz into the recursion relation Eq. ([3]), we 
obtain the following fixed-point condition: 



w{n) = w(m)w(n — 1 — m) 



(6) 



This equation can be solved by applying the discrete 
Laplace Transform w{z) = X)m=o w(m)z~ m . The so- 
lution z = 



, 4 — zw(0) has a branch-cut starting at 
z = l/4u>(0), with w(0) an undetermined constant. Af- 
ter applying an inverse Laplace Transform, one finds that 
w(n) actually has the same scaling form as the partition 
function of a molten globule: 



w(n) 



exp[-nln(l/4«j(0))] 



(7) 



However, the mathematical origin of the n~ 3 ^ 2 factor 
is here a combinatorial factor that reflects the different 
ways one can partition the open bubble between the two 
sub-trees. We may interpret £ ~ l/ln(l/4w(0)) as the 
characteristic size of a molten globule bubble at the root 
of the tree in the ordered phase. Numerical iteration 
of the recursion relations for W(k,n) for q = 150 and 
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q = 4 were found to be consistent with Eq. 0. If one 
uses W(N,n) « w{n)z^ in the remaining recursion re- 
lation Eq. with Eq. ©, one obtains the following 
self-consistency relation for the unknown w(0): 



w(0) 



2ni 



-G(z)w(l/z) 2 dz (8) 

z 



Here, G(z) is the discrete Laplace Transform of G(L) [Io| . 
which has a branch-cut that terminates at z = (1+2^/q) 2 . 
The integration contour in Eq. ([5]) must run inside an 
annulus in the complex plane that surrounds the origin 
passing the real axis outside the branch-cut of G{z) that 
terminates at (1 + 2^/q) 2 but inside the branch-cut of 
w(l/z) that starts at z = l/4w(0) . That means that 
the contour integral can only be carried out as long as 

The partition function develops a mathematical singular- 
ity when the two branch cuts merge, i.e., when Eq. ([9]) 
reduces to an equality. At that point, the partition 
w(n) ~ ^rcl/"^ °f the root bubble has the same form 
as the partition G(ri) « — for a molten globule 

of the same size. We can identify w(0) = a^qj^p as 
the stability limit of the groundstate. Note that the (low 
temperature) correlation length £ ~ l/ln(l/4w;(0)) can- 
not diverge at the stability limit. The critical value q c 
for q at the stability limited is now easily obtained by 
noting that w(0) is small compared to one. Expanding 
the argument of the contour integral in powers of w(0) 
leads to: 

(q^)- 1 ^w{0y r 2(l-iq)w(0) 2 +5(mq+2q 2 )w(0) 3 +. . . (10) 
If Eq. (|10p is combined with w(0) = 4n4 _l r y± one finds 



_ 4(1+2^)* 

that for q = 4, the singularity is at q sa 92.6. 

Surprisingly, the numerically computed free energy per 
site shown in Fig. 2 exhibits no singular dependence on q 
in that range. This is not inconsistent because w(n) only 
contributes a sub-leading term to the total free energy. 
On the other hand, the correlation length obtained from 
the pinching free energy appears to diverge near q c . We 
encountered however strong finite-size effects in the nu- 
merical solution of the recursion relations for q values in 
the range between 80 and 90 which make it difficult to 
numerically explore the critical properties in more detail. 
In addition, over that range of q values, our fixed-point 
scaling Ansatz appears not to be valid, at least for k val- 
ues less than 19. Instead, the reduced partition function 
scales as W(k, n)/W(k, 0) w N{k)g(n/N{k)) with g(x) a 
scaling function that is nearly linear for small values of x. 

In summary, a branched RNA molecule in the form 
of a Cayley Tree undergoes a phase transition from the 
branched groundstate to a molten globule phase if one 
reduces the energetic bias for the groundstate. The sta- 
bility of the branched groundstate against thermal fluctu- 



ations is significantly less than that of the linear ground- 
state. Branching does not produce mean-field critical be- 
havior but, instead, smears out the specific heat anomaly 
that characterizes systems with a linear groundstate. On 
the "high-temperature" side of the melting transition, 
numerical solution of the recursion relation produces a 
diverging correlation length. We showed - analytically - 
that on the low temperature side the designed ground- 
state becomes unstable at a critical point where the free 
energy develops a mathematical singularity not associ- 
ated with a divergence of the correlation length. 

Experimental studies comparing the melting charac- 
teristics of large, branched RNA molecules with that of 
linear, unbranched molecules that could probe this exotic 
form of melting have not yet been carried out but such 
systems would be fascinating laboratories for statistical 
mechanics. An important question in this respect would 
be the role of excluded volume interactions and of "ter- 
tiary" pairing interactions, i.e., pairing interactions that 
introduce, for example, pseudo-knots. Excluded volume 
interactions in general tend to suppress thermal fluctua- 
tions and possibly could restore the thermodynamic sin- 
gularity in the free energy per site that was encountered 
for linear molecules. Tertiary interaction could have the 
effect of turning a branched, secondary template into a 
three dimensional gel-like structure, in which case the 
transition to the molten-globule state could resemble the 
melting transition of a bulk solid material. 
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